package com.cyf.utils.sunpos;

public class SolPos {

    /*============================================================================
    *
    *     Define the function codes
    *
    *----------------------------------------------------------------------------*/
    public int L_DOY = 0x0001;
    public int L_GEOM = 0x0002;
    public int L_ZENETR = 0x0004;
    public int L_SSHA = 0x0008;
    public int L_SBCF = 0x0010;
    public int L_TST = 0x0020;
    public int L_SRSS = 0x0040;
    public int L_SOLAZM = 0x0080;
    public int L_REFRAC = 0x0100;
    public int L_AMASS = 0x0200;
    public int L_PRIME = 0x0400;
    public int L_TILT = 0x0800;
    public int L_ETR = 0x1000;
    public int L_ALL = 0xFFFF;

    /*============================================================================
*
*     Define the bit-wise masks for each function
*
*----------------------------------------------------------------------------*/
    int S_DOY = (L_DOY);
    int S_GEOM = (L_GEOM | S_DOY);
    int S_ZENETR = (L_ZENETR | S_GEOM);
    int S_SSHA = (L_SSHA | S_GEOM);
    int S_SBCF = (L_SBCF | S_SSHA);
    int S_TST = (L_TST | S_GEOM);
    int S_SRSS = (L_SRSS | S_SSHA | S_TST);
    int S_SOLAZM = (L_SOLAZM | S_ZENETR);
    int S_REFRAC = (L_REFRAC | S_ZENETR);
    int S_AMASS = (L_AMASS | S_REFRAC);
    int S_PRIME = (L_PRIME | S_AMASS);
    int S_TILT = (L_TILT | S_SOLAZM | S_REFRAC);
    int S_ETR = (L_ETR | S_REFRAC);
    int S_ALL = (L_ALL);

    /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    *
    * Structures defined for this module
    *
    *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    class trigdata /* used to pass calculated values locally */ {
        double cd;       /* cosine of the declination */
        double ch;       /* cosine of the hour angle */
        double cl;       /* cosine of the latitude */
        double sd;       /* sine of the declination */
        double sl;       /* sine of the latitude */
    }

    ;

    double degrad = 57.295779513; /* converts from radians to degrees */
    double raddeg = 0.0174532925; /* converts from degrees to radians */

    int[][] month_days = new int[][]{{0, 0, 31, 59, 90, 120, 151,
            181, 212, 243, 273, 304, 334},
            {0, 0, 31, 60, 91, 121, 152,
                    182, 213, 244, 274, 305, 335}};

    protected void geometry(PosData pdat) {
        double bottom;      /* denominator (bottom) of the fraction */
        double c2;          /* cosine of d2 */
        double cd;          /* cosine of the day angle or delination */
        double d2;          /* pdat->dayang times two */
        double delta;       /* difference between current year and 1949 */
        double s2;          /* sine of d2 */
        double sd;          /* sine of the day angle */
        double top;         /* numerator (top) of the fraction */
        int leap;        /* leap year counter */

  /* Day angle */
      /*  Iqbal, M.  1983.  An Introduction to Solar Radiation.
            Academic Press, NY., page 3 */
        pdat.dayang = 360.0 * (pdat.daynum - 1) / 365.0;

    /* Earth radius vector * solar constant = solar energy */
        /*  Spencer, J. W.  1971.  Fourier series representation of the
            position of the sun.  Search 2 (5), page 172 */
        sd = sin(raddeg * pdat.dayang);
        cd = cos(raddeg * pdat.dayang);
        d2 = 2.0 * pdat.dayang;
        c2 = cos(raddeg * d2);
        s2 = sin(raddeg * d2);

        pdat.erv = 1.000110 + 0.034221 * cd + 0.001280 * sd;
        pdat.erv += 0.000719 * c2 + 0.000077 * s2;

    /* Universal Coordinated (Greenwich standard) time */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.utime =
                pdat.hour * 3600.0 +
                        pdat.minute * 60.0 +
                        pdat.second -
                        (double) pdat.interval / 2.0;
        pdat.utime = pdat.utime / 3600.0 - pdat.timezone;

    /* Julian Day minus 2,400,000 days (to eliminate roundoff errors) */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */

    /* No adjustment for century non-leap years since this function is
       bounded by 1950 - 2050 */
        delta = pdat.year - 1949;
        leap = (int) (delta / 4.0);
        pdat.julday =
                32916.5 + delta * 365.0 + leap + pdat.daynum + pdat.utime / 24.0;

    /* Time used in the calculation of ecliptic coordinates */
    /* Noon 1 JAN 2000 = 2,400,000 + 51,545 days Julian Date */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.ectime = pdat.julday - 51545.0;

    /* Mean longitude */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.mnlong = 280.460 + 0.9856474 * pdat.ectime;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
        pdat.mnlong -= 360.0 * (int) (pdat.mnlong / 360.0);
        if (pdat.mnlong < 0.0)
            pdat.mnlong += 360.0;

    /* Mean anomaly */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.mnanom = 357.528 + 0.9856003 * pdat.ectime;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
        pdat.mnanom -= 360.0 * (int) (pdat.mnanom / 360.0);
        if (pdat.mnanom < 0.0)
            pdat.mnanom += 360.0;

    /* Ecliptic longitude */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.eclong = pdat.mnlong + 1.915 * sin(pdat.mnanom * raddeg) +
                0.020 * sin(2.0 * pdat.mnanom * raddeg);

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
        pdat.eclong -= 360.0 * (int) (pdat.eclong / 360.0);
        if (pdat.eclong < 0.0)
            pdat.eclong += 360.0;

    /* Obliquity of the ecliptic */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */

    /* 02 Feb 2001 SMW corrected sign in the following line */
/*  pdat.ecobli = 23.439 + 4.0e-07 * pdat.ectime;     */
        pdat.ecobli = 23.439 - 4.0e-07 * pdat.ectime;

    /* Declination */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.declin = degrad * asin(sin(pdat.ecobli * raddeg) *
                sin(pdat.eclong * raddeg));

    /* Right ascension */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        top = cos(raddeg * pdat.ecobli) * sin(raddeg * pdat.eclong);
        bottom = cos(raddeg * pdat.eclong);

        pdat.rascen = degrad * atan2(top, bottom);

    /* (make it a positive angle) */
        if (pdat.rascen < 0.0)
            pdat.rascen += 360.0;

    /* Greenwich mean sidereal time */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.gmst = 6.697375 + 0.0657098242 * pdat.ectime + pdat.utime;

    /* (dump the multiples of 24, so the answer is between 0 and 24) */
        pdat.gmst -= 24.0 * (int) (pdat.gmst / 24.0);
        if (pdat.gmst < 0.0)
            pdat.gmst += 24.0;

    /* Local mean sidereal time */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.lmst = pdat.gmst * 15.0 + pdat.longitude;

    /* (dump the multiples of 360, so the answer is between 0 and 360) */
        pdat.lmst -= 360.0 * (int) (pdat.lmst / 360.0);
        if (pdat.lmst < 0.)
            pdat.lmst += 360.0;

    /* Hour angle */
        /*  Michalsky, J.  1988.  The Astronomical Almanac's algorithm for
            approximate solar position (1950-2050).  Solar Energy 40 (3),
            pp. 227-235. */
        pdat.hrang = pdat.lmst - pdat.rascen;

    /* (force it between -180 and 180 degrees) */
        if (pdat.hrang < -180.0)
            pdat.hrang += 360.0;
        else if (pdat.hrang > 180.0)
            pdat.hrang -= 360.0;
    }

    void zen_no_ref(PosData pdat, Trigdata tdat) {
        double cz;          /* cosine of the solar zenith angle */

        localtrig(pdat, tdat);
        cz = tdat.sd * tdat.sl + tdat.cd * tdat.cl * tdat.ch;

    /* (watch out for the roundoff errors) */
        if (fabs(cz) > 1.0) {
            if (cz >= 0.0)
                cz = 1.0;
            else
                cz = -1.0;
        }

        pdat.zenetr = acos(cz) * degrad;

    /* (limit the degrees below the horizon to 9 [+90 -> 99]) */
        if (pdat.zenetr > 99.0)
            pdat.zenetr = 99.0;

        pdat.elevetr = 90.0 - pdat.zenetr;
    }

    public long S_solpos(PosData pdat) {
        Trigdata tdat = new Trigdata();

  /* initialize the trig structure */
        tdat.sd = -999.0; /* flag to force calculation of trig data */
        tdat.cd = 1.0;
        tdat.ch = 1.0; /* set the rest of these to something safe */
        tdat.cl = 1.0;
        tdat.sl = 1.0;

        doy2dom(pdat);                /* convert input doy to month-day */

        geometry(pdat);               /* do basic geometry calculations */

        zen_no_ref(pdat, tdat);

        sazm(pdat, tdat);
        refrac(pdat);

        return 0;
    }

    int LOCAL_TEMP = 25; //当地温度用来校正
    int LOCAL_PRESS = 1006;//当地压强用来校正（0~2000）

    void refrac(PosData pdat) {
        double prestemp;    /* temporary pressure/temperature correction */
        double refcor;      /* temporary refraction correction */
        double tanelev;     /* tangent of the solar elevation angle */

    /* If the sun is near zenith, the algorithm bombs; refraction near 0 */
        if (pdat.elevetr > 85.0)
            refcor = 0.0;

    /* Otherwise, we have refraction */
        else {
            tanelev = tan(raddeg * pdat.elevetr);
            if (pdat.elevetr >= 5.0)
                refcor = 58.1 / tanelev -
                        0.07 / (pow(tanelev, 3)) +
                        0.000086 / (pow(tanelev, 5));
            else if (pdat.elevetr >= -0.575)
                refcor = 1735.0 +
                        pdat.elevetr * (-518.2 + pdat.elevetr * (103.4 +
                                pdat.elevetr * (-12.79 + pdat.elevetr * 0.711)));
            else
                refcor = -20.774 / tanelev;

            prestemp =
                    (LOCAL_PRESS * 283.0) / (1013.0 * (273.0 + LOCAL_TEMP));
            refcor *= prestemp / 3600.0;
        }

    /* Refracted solar elevation angle */
        pdat.elevref = pdat.elevetr + refcor;

    /* (limit the degrees below the horizon to 9) */
        if (pdat.elevref < -9.0)
            pdat.elevref = -9.0;


    }


    void sazm(PosData pdat, Trigdata tdat) {
        double ca;          /* cosine of the solar azimuth angle */
        double ce;          /* cosine of the solar elevation */
        double cecl;        /* ( ce * cl ) */
        double se;          /* sine of the solar elevation */

        localtrig(pdat, tdat);
        ce = cos(raddeg * pdat.elevetr);
        se = sin(raddeg * pdat.elevetr);

        pdat.azim = 180.0;
        cecl = ce * tdat.cl;
        if (fabs(cecl) >= 0.001) {
            ca = (se * tdat.sl - tdat.sd) / cecl;
            if (ca > 1.0)
                ca = 1.0;
            else if (ca < -1.0)
                ca = -1.0;

            pdat.azim = 180.0 - acos(ca) * degrad;
            if (pdat.hrang > 0)
                pdat.azim = 360.0 - pdat.azim;
        }
    }

    /*============================================================================
*    Local Void function localtrig
*
*    Does trig on internal variable used by several functions
*----------------------------------------------------------------------------*/
    void localtrig(PosData pdat, Trigdata tdat) {

        if (tdat.sd < -900.0)  /* sd was initialized -999 as flag */ {
            tdat.sd = 1.0;  /* reflag as having completed calculations */
            tdat.cd = cos(raddeg * pdat.declin);

            tdat.ch = cos(raddeg * pdat.hrang);

            tdat.cl = cos(raddeg * pdat.latitude);

            tdat.sd = sin(raddeg * pdat.declin);

            tdat.sl = sin(raddeg * pdat.latitude);
        }
    }


    void doy2dom(PosData pdat) {
        int imon;  /* Month (month_days) array counter */
        int leap;  /* leap year switch */

    /* Set the leap year switch */
        if (((pdat.year % 4) == 0) &&
                (((pdat.year % 100) != 0) || ((pdat.year % 400) == 0)))
            leap = 1;
        else
            leap = 0;

    /* Find the month */
        imon = 12;
        while (pdat.daynum <= month_days[leap][imon])
            --imon;

    /* Set the month and day of month */
        pdat.month = imon;
        pdat.day = pdat.daynum - month_days[leap][imon];
    }

    private double sin(double v) {
        return Math.sin(v);
    }

    private double cos(double v) {
        return Math.cos(v);
    }

    private double acos(double v) {
        return Math.acos(v);
    }

    private double asin(double v) {
        return Math.asin(v);
    }

    private double atan2(double top, double bottom) {
        return Math.atan2(top, bottom);
    }

    private double fabs(double v) {
        return Math.abs(v);
    }

    private double tan(double v) {
        return Math.tan(v);
    }

    private double pow(double a, double b) {
        return Math.pow(a, b);
    }
}
